#clean up
rm(list=ls())

#Set Working Directory: REVISE THIS TO WHERE YOU SAVE THE DATA.
setwd('/Volumes/MONOGAN/mediaCount/data/')

#Load Libraries
library(MASS)
library(TSA)
library(foreign)
library(lattice)

#Load Patrick Brandt's code for a Poisson autoregressive model:
source('http://www.utdallas.edu/~pbrandt/code/pests.r')

########PEAKE & ESHBAUGH-SOHA 2008, TABLE 2##############
pes.2<-read.csv('PESenergy.csv')

pdf("pesHist.pdf",family='Times',width=6,height=6,pointsize=10)
#postscript("pesHist.eps",family='ComputerModern',width=6,height=6,pointsize=10)
plot(table(pes.2$Energy),type='h',axes=F,xlab="Energy (TV News)",ylab="Frequency",xlim=c(0,225))
axis(1,at=seq(0,230,by=10))
axis(2)
abline(h=0,col='gray60')
dev.off()

###Negative Binomial Regression###
pes.table.2.nb<-glm.nb(pes.2$Energy~pes.2$rmn1173+pes.2$grf0175+pes.2$grf575+pes.2$jec477+pes.2$jec1177+pes.2$jec479+pes.2$embargo+pes.2$hostages+pes.2$oilc+pes.2$Approval+pes.2$Unemploy); summary(pes.table.2.nb)

#correlation between true and fitted
cor(pes.2$Energy[-1],pes.table.2.nb$fitted.values[-1])
cor(pes.2$Energy[-1],pes.table.2.nb$fitted.values[-1])^2

#LB
Box.test(pes.table.2.nb$residuals,16,'Ljung-Box')

#RMSE
res<-pes.2$Energy-pes.table.2.nb$fitted.values
sqrt(sum(res[-1]^2))

#Jarque-Bera
jarque.bera.test(pes.table.2.nb$residuals)

#long-run effects
100*(exp(pes.table.2.nb$coefficients)-1)

###Self-Identify AR process and fit PAR model###
#standardized residuals from NB reg
pes.2.std<-(pes.2$Energy-pes.table.2.nb$fitted.values)/sqrt(pes.table.2.nb$fitted.values)

#ACF/PACF
acf(pes.2.std)
pacf(pes.2.std)

##Load New Code to Estimate Poisson Autoregressive Model with Bounds on Autoregressive Parameter##
#Note: The number of upper and lower bounds presented is model-sensitive.
Parp<-function(formula, p=1, parp.init=rep(0.1,p), init.param=NULL, ...)
  {
    model <- model.frame(formula)
    y <- model.extract(model, "response")
    x <- model.matrix(formula)
    k <- ncol(x)
    xnames <- dimnames(x)[[2]]
    colnames(x) <- xnames

    # Do some work to get a vector of starting values.
    if(is.null(init.param))
       {
         init.param <- glm(formula, family=poisson())$coefficients
       }

    # Estimate the PAR(p) model using L-BFGS-B
    fitted.param<-optim(c(parp.init,init.param), parpllf, method = "L-BFGS-B", hessian=TRUE,
#The next two lines create lower and upper bounds for constrained optimization.
#In effect, the first parameter (autoregressive) is all we restrict.
    						lower=c(0,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf,-Inf), 
    						upper=c(1,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf,Inf),
                        control=list(maxit=1000,
                          trace=2, fnscale=-(length(y))),
                        y=y,x=x,...)

    # Retreive the parameters and llf values
    param<-fitted.param$par
    llf <- fitted.param$value

    # Estimate the standard errors and compute z scores
    covar <- -solve(fitted.param$hessian)
    se<-sqrt(diag(covar))
    z<-param/se

    k <- ncol(x)
    k.plus.p <- length(se)
    num.p <- k.plus.p - k
    fitted.par.filter <- parpfilter(y,x,param[1:num.p],param[(num.p+1):k.plus.p])

    residuals <- y - fitted.par.filter[,"m"]
    # Print the results
    coefs <- t(matrix(rbind(param,se,z),nrow=3));
    rownames(coefs) <- c(rep("rho",num.p),colnames(x))
    colnames(coefs) <- c("Parameters","Std. Errors","Z-score")
    aic <- (-2*llf)+(2*(length(param)-1))
    dof <- length(y) - k.plus.p
    R <- diag(1, nrow=p, ncol=k.plus.p)
    wald.test <- wald(matrix(param, ncol=1), R, covar, length(y))
    wald.p <- wald.test$Wald.p
    wald.test <- wald.test$Wald.stat

    fit <- list(coefs=coefs,
                optim.param=fitted.param,
                param=param,
                covar=covar,
                std.err=se,
                z=z,
                residuals=residuals,
                llf=llf,
                aic=aic,
                dof=dof,
                wald.test=wald.test,
                wald.p=wald.p,
                k=k,
                p=num.p,
                y=y,
                x=x,
                call=match.call(),
                obj.type=c("PAR(p) regression output"))
    class(fit) <- c("pests.output")
    return(fit)
  }

##PAR model##
pes.table.2.par1<-Parp(pes.2$Energy~pes.2$rmn1173+pes.2$grf0175+pes.2$grf575+pes.2$jec477+pes.2$jec1177+pes.2$jec479+pes.2$embargo+pes.2$hostages+pes.2$oilc+pes.2$Approval+pes.2$Unemploy, p=1); summary(pes.table.2.par1)
pes.2.parp.std<-(pes.table.2.par1$residuals)/sqrt(pes.2$Energy-pes.table.2.par1$residuals)
acf(pes.2.parp.std,26)
pacf(pes.2.parp.std,26)
Box.test(pes.2.parp.std,26,'Ljung-Box')
Box.test(pes.2.parp.std,16,'Ljung-Box')


#correlation between true and fitted
par.fit<-pes.2$Energy-pes.table.2.par1$residuals
cor(pes.2$Energy[-1],par.fit[-1])
cor(pes.2$Energy[-1],par.fit[-1])^2

#LB
Box.test(pes.2.parp.std,16,'Ljung-Box')
Box.test(pes.table.2.par1$residuals,16,'Ljung-Box')

#RMSE
sqrt(sum(pes.table.2.par1$residuals[-1]^2))

#Jarque-Bera
jarque.bera.test(pes.table.2.par1$residuals)

#long-run effects
parp.multipliers(pes.table.2.par1,xvec=matrix(c(1,0,0,0,0,0,0,1,0,36.27,48.62,6.65),1,pes.table.2.par1$k))
100*((-0.083180747840287810879/par.fit[57]))
parp.multipliers(pes.table.2.par1,xvec=matrix(c(1,0,0,0,0,0,0,0,1,36.27,48.62,6.65),1,pes.table.2.par1$k))
100*((-0.72228783490700942949/par.fit[130]))
parp.multipliers(pes.table.2.par1,xvec=matrix(c(1,0,0,0,0,0,0,0,0,36.27,48.62,6.65),1,pes.table.2.par1$k))
100*((-4.4475243124047203480/34.32))
100*((-0.74222777669808015499/34.32))
100*((-2.35074961571743212829/34.32))

###Transfer Function###
pes.table.2<-arimax(pes.2$Energy, order=c(1,0,0), xtransf=cbind(pes.2$rmn1173,pes.2$grf0175,pes.2$grf575,pes.2$jec477,pes.2$jec1177,pes.2$jec479,pes.2$embargo,pes.2$hostages,pes.2$oilc,pes.2$Approval,pes.2$Unemploy), transfer=list(c(1,0),c(0,0),c(0,0),c(0,0),c(0,0),c(0,0),c(0,0),c(0,0),c(0,0),c(0,0),c(0,0))); pes.table.2
acf(pes.table.2$residuals)
pacf(pes.table.2$residuals)
Box.test(pes.table.2$residuals,16,'Ljung-Box')

#correlation between true and fitted
tf.fit<-pes.2$Energy-pes.table.2$residuals
cor(pes.2$Energy[-1],tf.fit[-1])
cor(pes.2$Energy[-1],tf.fit[-1])^2

#LB
Box.test(pes.table.2$residuals,16,"Ljung-Box")

#RMSE
sqrt(sum(pes.table.2$residuals[-1]^2))

#Jarque-Bera
jarque.bera.test(pes.table.2$residuals)

#long-run effects
100*((26.579516926218094142/tf.fit[57]))
100*((-15.515078005996985056/tf.fit[130]))
100*((0.27360887233532704688/34.32))
100*((-0.27692805557702265817/34.32))
100*((0.081216765120832964953/34.32))

###OLS with 1 lag. (Koyck model.)###
t.Energy<-ts(pes.2$Energy)
t.rmn1173<-ts(pes.2$rmn1173)
t.grf0175<-ts(pes.2$grf0175)
t.grf575<-ts(pes.2$grf575)
t.jec477<-ts(pes.2$jec477)
t.jec1177<-ts(pes.2$jec1177)
t.jec479<-ts(pes.2$jec479)
t.embargo<-ts(pes.2$embargo)
t.hostages<-ts(pes.2$hostages)
t.oilc<-ts(pes.2$oilc)
t.Approval<-ts(pes.2$Approval)
t.Unemploy<-ts(pes.2$Unemploy)
l.Energy<-lag(t.Energy,-1)
pes.2.b<-ts.union(t.Energy,l.Energy,t.rmn1173,t.grf0175,t.grf575,t.jec477,t.jec1177,t.jec479,t.embargo,t.hostages,t.oilc,t.Approval,t.Unemploy)
pes.table.2.koyck<-lm(t.Energy~l.Energy+t.rmn1173+t.grf0175+t.grf575+t.jec477+t.jec1177+t.jec479+t.embargo+t.hostages+t.oilc+t.Approval+t.Unemploy, data=pes.2.b); summary(pes.table.2.koyck)

#correlation between true and fitted
cor(pes.2$Energy[-1],pes.table.2.koyck$fitted.values)
cor(pes.2$Energy[-1],pes.table.2.koyck$fitted.values)^2

#LB
Box.test(pes.table.2.koyck$residuals,16,'Ljung-Box')

#RMSE
sqrt(sum(pes.table.2.koyck$residuals^2))

#Jarque-Bera
jarque.bera.test(pes.table.2.koyck$residuals)

#long-run effects
10.540609533265392628/(1-0.739229585491644747)
100*((40.421033011502402132/pes.table.2.koyck$fitted.values[57]))
-2.514122100581635877/(1-0.739229585491644747)
100*((-9.6411324318429603153/pes.table.2.koyck$fitted.values[57]))
-1.141710045498506076/(1-0.739229585491644747)
100*((-4.3782192379876931909/34.32))
-0.154375415799812010/(1-0.739229585491644747)
100*((-0.5919974322657132193/34.32))
-0.886553166635491885/(1-0.739229585491644747)
100*((-3.3997459731272012817/34.32))

###Log-Normal Model###
log.Energy<-ts(log(pes.2$Energy+.1))
l.log.Energy<-lag(log.Energy,-1)
pes.3.b<-ts.union(log.Energy,l.log.Energy,t.rmn1173,t.grf0175,t.grf575,t.jec477,t.jec1177,t.jec479,t.embargo,t.hostages,t.oilc,t.Approval,t.Unemploy)
pes.table.2.lognorm<-lm(log.Energy~l.log.Energy+t.rmn1173+t.grf0175+t.grf575+t.jec477+t.jec1177+t.jec479+t.embargo+t.hostages+t.oilc+t.Approval+t.Unemploy, data=pes.3.b); summary(pes.table.2.lognorm)

#correlation between true and fitted
cor(pes.2$Energy[-1],exp(pes.table.2.lognorm$fitted.values))
cor(pes.2$Energy[-1],exp(pes.table.2.lognorm$fitted.values))^2

#LB
Box.test(pes.table.2.lognorm$residuals,16,'Ljung-Box')

#RMSE
res<-pes.2$Energy[-1]-exp(pes.table.2.lognorm$fitted.values)
sqrt(sum(res^2))

#Jarque-Bera
jarque.bera.test(pes.table.2.lognorm$residuals)

#long-run effects
100*(exp(0.4720124895621945482/(1-0.5513763598916749942))-1)
100*(exp(0.0269420599912406702/(1-0.5513763598916749942))-1)
100*(exp(-0.1189758928677158711/(1-0.5513763598916749942))-1)
100*(exp(-0.0143443965953866624/(1-0.5513763598916749942))-1)
100*(exp(-0.0319464425779800720/(1-0.5513763598916749942))-1)


########DYNAMIC FIGURE OF THE EFFECT OF NIXON'S 11/73 SPEECH########
###Clarify-like simulations###
pes.2.sim.par<-mvrnorm(n=1000,mu=pes.table.2.par1$coefs[,1],Sigma=pes.table.2.par1$covar,empirical=TRUE)

pes.2.results.nb<-summary(pes.table.2.nb)
pes.2.sim.nb<-mvrnorm(n=1000,mu=pes.2.results.nb$coefficients[,1],Sigma=pes.2.results.nb$cov.unscaled,empirical=TRUE)

pes.2.results.koyck<-summary(pes.table.2.koyck)
pes.2.covar.koyck<-pes.2.results.koyck$cov.unscaled*(pes.2.results.koyck$sigma^2)
pes.2.sim.koyck<-mvrnorm(n=1000,mu=pes.2.results.koyck$coefficients[,1],Sigma=pes.2.covar.koyck,empirical=TRUE)

pes.2.results.lognorm<-summary(pes.table.2.lognorm)
pes.2.covar.lognorm<-pes.2.results.lognorm$cov.unscaled*(pes.2.results.lognorm$sigma^2)
pes.2.sim.lognorm<-mvrnorm(n=1000,mu=pes.2.results.lognorm$coefficients[,1],Sigma=pes.2.covar.lognorm,empirical=TRUE)

pes.2.sim.tf<-mvrnorm(n=1000,mu=pes.table.2$coef,Sigma=pes.table.2$var.coef,empirical=TRUE)

###Transfer function###
#Start with predictions from reported coefficients.
coef<-pes.table.2$coef[c(4:14,2)]
#c(162.8547,30.2612,-8.7277,29.8458,-6.7220,-18.5718,26.5795,-15.5151,.2736,-.2769,.0812,35.2992)
nixon.initial<-pes.table.2$coef[4]
#162.8547
nixon.decay<-pes.table.2$coef[3]
#.6167

#Nixon onsets at wave 1. Waves -1-10
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1)
x1<-c(1,0,0,0,0,0,1,0,37.1,30,4.6,1)
y00<-coef%*%x00
y0<-coef%*%x0
y1<-coef%*%x1
y2<-coef%*%x2+nixon.decay*nixon.initial
y3<-coef%*%x3+(nixon.decay^2)*nixon.initial
y4<-coef%*%x4+(nixon.decay^3)*nixon.initial
y5<-coef%*%x5+(nixon.decay^4)*nixon.initial
y6<-coef%*%x6+(nixon.decay^5)*nixon.initial
y7<-coef%*%x7+(nixon.decay^6)*nixon.initial
y8<-coef%*%x8+(nixon.decay^7)*nixon.initial
y9<-coef%*%x9+(nixon.decay^8)*nixon.initial
y10<-coef%*%x10+(nixon.decay^9)*nixon.initial

#Plot Forecast
Time<-c(-1:10)
mle.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=mle.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
coef<-pes.2.sim.tf[i,c(4:14,2)]
nixon.initial<-pes.2.sim.tf[i,4]
nixon.decay<-pes.2.sim.tf[i,3]
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1)
x1<-c(1,0,0,0,0,0,1,0,37.1,30,4.6,1)
forecast[i,1]<-y00<-coef%*%x00
forecast[i,2]<-y0<-coef%*%x0
forecast[i,3]<-y1<-coef%*%x1
forecast[i,4]<-y2<-coef%*%x2+nixon.decay*nixon.initial
forecast[i,5]<-y3<-coef%*%x3+(nixon.decay^2)*nixon.initial
forecast[i,6]<-y4<-coef%*%x4+(nixon.decay^3)*nixon.initial
forecast[i,7]<-y5<-coef%*%x5+(nixon.decay^4)*nixon.initial
forecast[i,8]<-y6<-coef%*%x6+(nixon.decay^5)*nixon.initial
forecast[i,9]<-y7<-coef%*%x7+(nixon.decay^6)*nixon.initial
forecast[i,10]<-y8<-coef%*%x8+(nixon.decay^7)*nixon.initial
forecast[i,11]<-y9<-coef%*%x9+(nixon.decay^8)*nixon.initial
forecast[i,12]<-y10<-coef%*%x10+(nixon.decay^9)*nixon.initial
}

#Lower and Upper Bounds
mle.p.5<-apply(forecast, 2, quantile,.05)
mle.p.95<-apply(forecast, 2, quantile,.95)

###PAR model###
#Start with predictions from reported coefficients.
rho<-pes.table.2.par1$coefs[1,1]
weight<-1-sum(rho)
beta<-pes.table.2.par1$coef[c(3:13,2),1]
#c(2.9099,1.1855,.7011,2.1101,-.4189,1.338,-.0129,-.0295,-.1766,-.0295,-.0933,11.6844)

#Nixon onsets at wave 1. Waves -1-10
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1)
x1<-c(1,0,0,0,0,0,1,0,37.1,30,4.6,1)
lag00<-c(44)
y00<-weight*exp(beta%*%x00)+rho%*%lag00
lag0<-c(44)
y0<-weight*exp(beta%*%x0)+rho%*%lag0
lag1<-c(44)
y1<-weight*exp(beta%*%x1)+rho%*%lag1
lag2<-c(y1)
y2<-weight*exp(beta%*%x2)+rho%*%lag2
lag3<-c(y2)
y3<-weight*exp(beta%*%x3)+rho%*%lag3
lag4<-c(y3)
y4<-weight*exp(beta%*%x4)+rho%*%lag4
lag5<-c(y4)
y5<-weight*exp(beta%*%x5)+rho%*%lag5
lag6<-c(y5)
y6<-weight*exp(beta%*%x6)+rho%*%lag6
lag7<-c(y6)
y7<-weight*exp(beta%*%x7)+rho%*%lag7
lag8<-c(y7)
y8<-weight*exp(beta%*%x8)+rho%*%lag8
lag9<-c(y8)
y9<-weight*exp(beta%*%x9)+rho%*%lag9
lag10<-c(y9)
y10<-weight*exp(beta%*%x10)+rho%*%lag10

#Plot Forecast
Time<-c(-1:10)
par.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=par.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
beta<-pes.2.sim.par[i,c(3:13,2)]
rho<-pes.2.sim.par[i,1]
weight<-1-sum(rho)
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1)
x1<-c(1,0,0,0,0,0,1,0,37.1,30,4.6,1)
lag00<-c(44)
forecast[i,1]<-y00<-weight*exp(beta%*%x00)+rho%*%lag00
lag0<-c(44)
forecast[i,2]<-y0<-weight*exp(beta%*%x0)+rho%*%lag0
lag1<-c(44)
forecast[i,3]<-y1<-weight*exp(beta%*%x1)+rho%*%lag1
lag2<-c(y1)
forecast[i,4]<-y2<-weight*exp(beta%*%x2)+rho%*%lag2
lag3<-c(y2)
forecast[i,5]<-y3<-weight*exp(beta%*%x3)+rho%*%lag3
lag4<-c(y3)
forecast[i,6]<-y4<-weight*exp(beta%*%x4)+rho%*%lag4
lag5<-c(y4)
forecast[i,7]<-y5<-weight*exp(beta%*%x5)+rho%*%lag5
lag6<-c(y5)
forecast[i,8]<-y6<-weight*exp(beta%*%x6)+rho%*%lag6
lag7<-c(y6)
forecast[i,9]<-y7<-weight*exp(beta%*%x7)+rho%*%lag7
lag8<-c(y7)
forecast[i,10]<-y8<-weight*exp(beta%*%x8)+rho%*%lag8
lag9<-c(y8)
forecast[i,11]<-y9<-weight*exp(beta%*%x9)+rho%*%lag9
lag10<-c(y9)
forecast[i,12]<-y10<-weight*exp(beta%*%x10)+rho%*%lag10
}

#Lower and Upper Bounds
par.p.5<-apply(forecast, 2, quantile,.05)
par.p.95<-apply(forecast, 2, quantile,.95)


###Negative binomial regression###
#Start with predictions from reported coefficients.
coef<-pes.2.results.nb$coefficients[c(2:12,1),1]
#c(.7223,.2882,-.2276,.966,.5732,1.1415,1.1409,.0894,-.2766,-.0321,-.077,15.2993)

#Nixon onsets at wave 1. Waves -1-10
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1)
x1<-c(1,0,0,0,0,0,1,0,37.1,30,4.6,1)
y00<-exp(coef%*%x00)
y0<-exp(coef%*%x0)
y1<-exp(coef%*%x1)
y2<-exp(coef%*%x2)
y3<-exp(coef%*%x3)
y4<-exp(coef%*%x4)
y5<-exp(coef%*%x5)
y6<-exp(coef%*%x6)
y7<-exp(coef%*%x7)
y8<-exp(coef%*%x8)
y9<-exp(coef%*%x9)
y10<-exp(coef%*%x10)

#Plot Forecast
Time<-c(-1:10)
nb.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=nb.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
coef<-pes.2.sim.nb[i,c(2:12,1)]
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1)
x1<-c(1,0,0,0,0,0,1,0,37.1,30,4.6,1)
forecast[i,1]<-y00<-exp(coef%*%x00)
forecast[i,2]<-y0<-exp(coef%*%x0)
forecast[i,3]<-y1<-exp(coef%*%x1)
forecast[i,4]<-y2<-exp(coef%*%x2)
forecast[i,5]<-y3<-exp(coef%*%x3)
forecast[i,6]<-y4<-exp(coef%*%x4)
forecast[i,7]<-y5<-exp(coef%*%x5)
forecast[i,8]<-y6<-exp(coef%*%x6)
forecast[i,9]<-y7<-exp(coef%*%x7)
forecast[i,10]<-y8<-exp(coef%*%x8)
forecast[i,11]<-y9<-exp(coef%*%x9)
forecast[i,12]<-y10<-exp(coef%*%x10)
}

#Lower and Upper Bounds
nb.p.5<-apply(forecast, 2, quantile,.05)
nb.p.95<-apply(forecast, 2, quantile,.95)


###Log-Normal Model###
coef<-pes.table.2.lognorm$coefficients[-2]
rho1<-pes.table.2.lognorm$coefficients[2]
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0,0,0,0,0,1,0,37.1,30,4.6)
x1<-c(1,1,0,0,0,0,0,1,0,37.1,30,4.6)
y00<-exp(coef%*%x00)*(44^rho1)
y0<-exp(coef%*%x0)*(44^rho1)
y1<-exp(coef%*%x1)*(44^rho1)
y2<-exp(coef%*%x2)*(y1^rho1)
y3<-exp(coef%*%x3)*(y2^rho1)
y4<-exp(coef%*%x4)*(y3^rho1)
y5<-exp(coef%*%x5)*(y4^rho1)
y6<-exp(coef%*%x6)*(y5^rho1)
y7<-exp(coef%*%x7)*(y6^rho1)
y8<-exp(coef%*%x8)*(y7^rho1)
y9<-exp(coef%*%x9)*(y8^rho1)
y10<-exp(coef%*%x10)*(y9^rho1)

#Plot Forecast
Time<-c(-1:10)
lognorm.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=lognorm.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
coef<-pes.2.sim.lognorm[i,-2]
rho1<-pes.2.sim.lognorm[i,2]
x00<-x0<-x2<-x3<-x4<-x5<-x6<-x7<-x8<-x9<-x10<-c(1,0,0,0,0,0,0,1,0,37.1,30,4.6)
x1<-c(1,1,0,0,0,0,0,1,0,37.1,30,4.6)
forecast[i,1]<-y00<-exp(coef%*%x00)*(44^rho1)
forecast[i,2]<-y0<-exp(coef%*%x0)*(44^rho1)
forecast[i,3]<-y1<-exp(coef%*%x1)*(44^rho1)
forecast[i,4]<-y2<-exp(coef%*%x2)*(y1^rho1)
forecast[i,5]<-y3<-exp(coef%*%x3)*(y2^rho1)
forecast[i,6]<-y4<-exp(coef%*%x4)*(y3^rho1)
forecast[i,7]<-y5<-exp(coef%*%x5)*(y4^rho1)
forecast[i,8]<-y6<-exp(coef%*%x6)*(y5^rho1)
forecast[i,9]<-y7<-exp(coef%*%x7)*(y6^rho1)
forecast[i,10]<-y8<-exp(coef%*%x8)*(y7^rho1)
forecast[i,11]<-y9<-exp(coef%*%x9)*(y8^rho1)
forecast[i,12]<-y10<-exp(coef%*%x10)*(y9^rho1)
}

#Lower and Upper Bounds
lognorm.p.5<-apply(forecast, 2, quantile,.05)
lognorm.p.95<-apply(forecast, 2, quantile,.95)


###OLS with 1 lag. (Koyck.)###
#Start with predictions from reported coefficients.
coef<-pes.2.results.koyck$coefficients[c(3:13,1:2),1]

#Nixon onsets at wave 1. Waves -1-10
x00<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,44)
y00<-coef%*%x00
x0<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,44)
y0<-coef%*%x0
x1<-c(1,0,0,0,0,0,1,0,37.1,30,4.6,1,44)
y1<-coef%*%x1
x2<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y1)
y2<-coef%*%x2
x3<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y2)
y3<-coef%*%x3
x4<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y3)
y4<-coef%*%x4
x5<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y4)
y5<-coef%*%x5
x6<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y5)
y6<-coef%*%x6
x7<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y6)
y7<-coef%*%x7
x8<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y7)
y8<-coef%*%x8
x9<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y8)
y9<-coef%*%x9
x10<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y9)
y10<-coef%*%x10

#Plot Forecast
Time<-c(-1:10)
ols.Forecast<-c(y00,y0,y1,y2,y3,y4,y5,y6,y7,y8,y9,y10)
plot(y=ols.Forecast,x=Time,type='o')

#Clarify-like confidence intervals
forecast<-matrix(NA,nrow=1000,ncol=12)

for (i in 1:1000){
coef<-pes.2.sim.koyck[i,c(3:13,1:2)]
x00<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,44)
forecast[i,1]<-y00<-coef%*%x00
x0<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,44)
forecast[i,2]<-y0<-coef%*%x0
x1<-c(1,0,0,0,0,0,1,0,37.1,30,4.6,1,44)
forecast[i,3]<-y1<-coef%*%x1
x2<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y1)
forecast[i,4]<-y2<-coef%*%x2
x3<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y2)
forecast[i,5]<-y3<-coef%*%x3
x4<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y3)
forecast[i,6]<-y4<-coef%*%x4
x5<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y4)
forecast[i,7]<-y5<-coef%*%x5
x6<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y5)
forecast[i,8]<-y6<-coef%*%x6
x7<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y6)
forecast[i,9]<-y7<-coef%*%x7
x8<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y7)
forecast[i,10]<-y8<-coef%*%x8
x9<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y8)
forecast[i,11]<-y9<-coef%*%x9
x10<-c(0,0,0,0,0,0,1,0,37.1,30,4.6,1,y9)
forecast[i,12]<-y10<-coef%*%x10
}

#Lower and Upper Bounds
ols.p.5<-apply(forecast, 2, quantile,.05)
ols.p.95<-apply(forecast, 2, quantile,.95)

###Plot All Forecasts Against Real Data###
pdf("pesTable2.pdf",family='Times',width=6,height=6,pointsize=10)
#postscript("pesTable2.eps",family='ComputerModern',width=6,height=6,pointsize=10)
Time<-c(-1:10)

plot(y=ols.Forecast,x=Time-.26,pch=23,ylim=c(0,845),xlim=c(-1,11),xlab="Time",ylab="Predicted Count",col='blue', axes=F)
arrows(x0=Time-.26,y0=ols.p.5,x1=Time-.26,y1=ols.p.95, code=3, angle=90, length=.0,lty=2,col='blue')
axis(2)
axis(1,at=c(0,2,4,6,8,10),c('Oct 73','Dec 73','Feb 74','Apr 74','June 74','Aug 74'))
box()

points(y=nb.Forecast,x=Time-.13,pch=20,col='red')
arrows(x0=Time-.13,y0=nb.p.5,x1=Time-.13,y1=nb.p.95, code=3, angle=90, length=.0,lty=3,col='red')

points(y=lognorm.Forecast,x=Time,pch=15,col='purple')
arrows(x0=Time,y0=lognorm.p.5,x1=Time,y1=lognorm.p.95, code=3, angle=90, length=.0,lty=5,col='purple')

points(y=mle.Forecast,x=Time+.13,pch=22,col='forestgreen')
arrows(x0=Time+.13,y0=mle.p.5,x1=Time+.13,y1=mle.p.95, code=3, angle=90, length=.0,lty=4,col='forestgreen')

points(y=par.Forecast,x=Time+.26,pch=24)
arrows(x0=Time+.26,y0=par.p.5,x1=Time+.26,y1=par.p.95, code=3, angle=90, length=.0,lty=1)

legend(x=6,y=800,legend=c('Koyck','Negative Binomial','Log-normal','Transfer Function','PAR'), lty=c(2,3,5,4,1), pch=c(23,20,15,22,24), col=c('blue','red','purple','forestgreen','black'))

abline(h=0, col='gray60')
dev.off()

